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Abstract Effective conductivity of a 2D composite corresponding to the regular 
hexagonal arrangement of superconducting disks is expressed in the form of a long 
series in the volume fraction of ideally conducting disks. According to our calcu¬ 
lations based on various re-summation techniques, both the threshold and critical 
index are obtained in good agreement with expected values. The critical amplitude 
is in the interval (5.14,5.24) that is close to the theoretical estimation 5.18. The 
next order (constant) term in the high concentration regime is calculated for the first 
time, and the best estimate is equal to (—6.229). Final formula is derived for the 
effective conductivity for arbitrary volume fraction. 
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1 Introduction 


We consider a two-dimensional composite corresponding to the regular hexagonal 
lattice arrangement of ideally conducting (superconducting) cylinders of radius a 
embedded into the matrix of a conducting material. The studies of the effective con¬ 
ductivity cr(x) of regular composites were pioneered by Maxwell lf36l l and Rayleigh 
03) . The results of these fundamental research remained limited to the lowest orders 
in x. Their work was extended in 02), resulting in rather good numerical solutions 
valid in much broader concentration intervals. 

The effective conductivity c(x) is an analytic function in x. In general case of a 
two-phase composite the so-called contrast parameter should be also included into 
consideration explicitly, see e.g. uni . We are interested in the case of a high contrast 
regular composites, when the conductivity of the inclusions is much larger than the 
conductivity of the host. I.e., the highly conducting inclusions are replaced by the 
ideally conducting inclusions with infinite conductivity. In this case the contrast 
parameter is equal to unity and remains implicit. The conductivity of the matrix is 
normalized by unity as well. From the phase interchange theorem lf33l it follows 
that in two-dimensions the superconductivity problem is dual to the conductivity 
problem, and the superconductivity critical index is equal to the conductivity index. 

Our study is restricted to the two-dimensional case which is still interesting, both 
for practical HE) and physical reasons [ 42l|47l l. Composite materials often consist 
of a uniform background-host reinforced by a large number (high concentration) of 
unidirectional rod -or fiber-like inclusions with high conductivity 0). 

On the other hand, two-dimensional regular hexagonal- arrayed composites 0), 
much closer resemble the two-dimensional random composites, than their respec¬ 
tive 3D counterparts do 03 - The tendency to order in the two-dimensional random 
system of disks, is a crucial feature in the theory of composites at high concentra¬ 
tions. 

Most strikingly it appears that the maximum volume fraction of -4= « 0.9069, 
is attained both for the regular hexagonal array of disks and for random (irregular) 
2D composites 02). 

A numerical study of the 2D hexagonal case can be found in 02). Their final 
formula (|T|) 

' ().075422.t 6 , i ' ^ 

1 - 1 . 06028 a 12 + 

compares rather well with numerical data of |42). Note that (JT]) diverges with critical 
exponent s = 1, as x —> 0.922351. This property on one hand makes the formula 
more accurate in the vicinity of a true critical point but, on the other hand, makes 
any comparison in the critical region meaningless. It remains rather accurate till 
x = 0.85, where the error is 0.47%. Forx = 0.905 the error is 52%. Expression 0 
was derived using only terms up to the 12th order in concentration. The expansion 
of 0 is characterized by a rather regular behavior of the coefficients. 
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O reg (x) = 1 + 2x + 2x 2 + 2x 3 + 2x 4 + 2x 5 + 2x 6 + 

2.15084x 7 + 2.30169x 8 + 2.45253x 9 + 2.60338x 10 + (2) 

2.75422x n +2.90506x 12 + <9(x 13 ). 

One can, in principle, collect the higher-order terms as well. However, such 
derivation of an additional terms can not be considered as consistent since it relies 
on the agreement with numerical results. It turns out though, that ([2} compares well 
with our results shown below, see ©. Except an immediate vicinity of the critical 
point, analytic-numeric approach of 0, is in a good agreement with the numerical 
results of ||42| . 

In a different limit of high concentrations Keller ll32l suggested a constructive 
asymptotic method for regular lattices, leading to very transparent, inverse square- 
root formula for the square array ll32l . Berlyand and Novikov (HD extended Keller’s 
method to the hexagonal array, 



Thus the critical amplitude A (pre-factor), is equal to A w 5.18. 

We will examine below this result for the critical amplitude from the perspec¬ 
tive of re-summation techniques suggested before for square regular arrays IH.By 
analogy with square lattice lf 38 l , we expect a constant correction in the asymptotic 
regime, 


cr ~ 


^3tt 3 / 2 


V2 


+ B, 


vT5 ' 


(4) 


where the correction term B can not be found in the literature, to the best of our 
knowledge. It will be calculated below by different methods. 

With account for such correction, the final universal formula valid for all possible 


concentrations from 0 to 




has the form 


<7 (x) = a{x) 


Pi*) 

2M’ 


(5) 


where 


a(x) = 


36.1415 


15.9909 


vT2 


— X 


Vl2 


— x —45.685 + 2.46148x, 


P(x) = (0.939152 + x) (1.38894 —2.16685x + x 2 )x 
(2.55367-0.836613x + x 2 )(2.08347+ 2.12786x+x 2 ) 


and 


Q(x) = (1.01215 +x)(1.61369 — 2 . 31669 x + x 2 )x 
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(6.51762 — 0.173965x+x 2 ) (4.88614 + 3.28716 x + x 2 ). 

The rest of the paper is organized as follows: in Section[2]we describe the essen¬ 
tials of the long series derivation. Section [3] applies various methods to the critical 
point calculation and compares the obtained results. In Section [4] the critical index 
and amplitude A are calculated. Section 0 where the most accurate formula for all 
volume fractions is derived, comparing the obtained predictions to numerical data. 
The amplitude B is calculated. Section[6]is concerned with interpolation with Pade 
approximants. Section [7] returns to discussion of the ansatz for construction of the 
starting approximation. Section Ogives unified approach to the square and hexago¬ 
nal lattices. Section [9] considers Dirichlet summation to extract the asymptotic be¬ 
havior of series coefficients. Section[10]derives the asymptotic formula by use of the 
lubrication theory. Section|TT]considers random composites related to the hexagonal 
lattice. Finally, Section[l2]concludes with a discussion of obtained results. 


2 Series for Hexagonal array of superconducting cylinders 


We proceed to the case of a hexagonal lattice of inclusions, where rather long expan¬ 
sions in concentration will be presented an analyzed systematically. The coefficients 
a n in the expansion of <r(x) = 1 + a nX n , are expressed through elliptic functions 
by exact formulas from lf39ll40B . Below, this expansion is presented in the truncated 
numerical form. 


a (x) = 1 + 2x + 2x 2 + 2x 3 + 2x 4 + 2x 5 + 2x 6 

+ 2.1508443464271876x 7 + 2.301688692854377x 8 
+ 2.452533039281566x 9 + 2.6033773857087543x 10 
+ 2.754221732135944x* 1 + 2.9050660785631326x 12 
+ 3.0674404324522926x 13 +3.2411917947659736‘x 14 
+ 3.426320165504177x ls + 3.6228255446669055x 16 
+ 3.8307079322541555x n + 4.049967328265928x 18 
+ 4.441422739726373x 19 + 4.845994396051242x 20 
+ 5.264540375940583x 21 + 5.69791875809444X 22 
+ 6.146987621212864X 23 + 6.6126050439959X 24 
+ 7.135044602470776x 25 + 7.700073986554016x 26 
+ 0(x 27 ). (6) 

The first twelve coefficients of © and the Taylor expansions of (©) coincide. The 
next coefficients can be calculated by exact formulas from I39ll40) . This requires 
use of the double precision and perhaps a power computer, not a standard laptop. 

Since we are dealing with the limiting case of a perfectly conducting inclusions 
when the conductivity of inclusions tends to infinity, the effective conductivity is 
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also expected to tend to infinity as a power-law, as the concentration jt tends to the 
maximal value x c for the hexagonal array, 

ff(x) ~ A(x c — + (7) 

The critical superconductivity index (exponent) s is believed to be 1 /2 for all lattices 
IH. For sake of exploring how consistent are various resummation techniques, we 
will calculate the index. The critical amplitudes A and B are unknown non-universal 
parameters to be calculated below as well. 

The problem of interest can be formulated mathematically as follows. Given the 
polynomial approximation (( 6 ji of the function cr(x), to estimate the convergence 
radius x c of the Taylor series <r(x); to determine critical index s and amplitudes A. B 
of the asymptotically equivalent approximation ([7} near x = x c . 

When such extrapolation problem is solved, we proceed to solve an interpolation 
problem of matching the two asymptotic expressions for the conductivity and derive 
interpolation formula for all concentrations. 


3 Critical Point 
3.1 Pade approximants 

Probably the simplest and direct way to extrapolate, is to apply the Pade approx¬ 
imants P n ,m{x), which is nothing else but ratio of the two polynomials P„(x) and 
P m (x) of the order n and m, respectively. The coefficients are derived directly from 
the coefficients of the given power series |6||44] from the requirement of asymptotic 
equivalence to the given series or function /(x). When there is a need to stress the 
last point, we simply write PadeApproximant[f[x\,n,m\. 

In order to estimate the position of a critical point, let us apply the diagonal Pade 
approximants. 


Pi,i{x) 


mix+ 1 
ii\x T 1 ’ 


Pi,i{x) 


ni 2 X 2 + m\x + 1 
112 X 2 + n ix + 1 


( 8 ) 


Pade approximants locally are the best rational approximations of power series. 
Their poles determine singular points of the approximated functions fi44l l6l. Cal¬ 
culations with Pade approximants are straightforward and can be performed with 
Mathematical. They do not require any preliminary knowledge of the critical in¬ 
dex and we have to find the position of a simple pole. In the theory of periodic 
2D composites H]|25)|45]|, their application is justifiable rigorously away from the 
square-root singularity and from the high-contrast limit. 

There is a convergence within the approximations for the critical point generated 
by the sequence of Pade approximants, corresponding to their order increasing: 
xi = 1, x '2 = 1, X 3 = 1, X 4 — n.a„ X 5 — n.a., X(, = 0.945958, X 7 = 0.945929, xs = 
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0.947703, xg = 0.946772, x 10 = 0.942378, x n = 0.945929, x 12 = 0.945959, x n = 
0.920878. 

The main body of the approximations is well off the exact value. The percentage 
error given by the last/best approximant in the sequence equals to 1.5413%. If only 
the first row of the Pade table is studied El, then the best estimate is equal to 
0.929867, close to the estimates with the diagonal sequence. 

We suggest that further increase in accuracy is limited by triviality, or “flatness” 
of the coefficients values in six starting orders of @. Consider another sequence 
of approximants, when diagonal Pade approximants are multiplied with Clausius- 
Mossotti-type expression, 


pH \ = d 1 -*) (1+'»!-*•) , 
1 (1 +x) (1 +n\x) ’ 

(1 —x) (1 + m\x + m 2 x 2 ) 

( 1 +x) (l+nix+n 2 x 2 ) ’ 


(9) 


The transformation which lifts the flatness, does improve convergence of the se¬ 
quence of approximations for the threshold, 

x~i = 0.94568, xg — n.a., xg = 0.948299, xio = 0.9287, xn = 0.945681, x\ 2 = 
0.89793, X 13 = 0.903517. The percentage error given by the last approximant in 
the sequence equals —0.373%. 

In order to judge the quality of the latter estimate, let us try highly recommended 
D — Log Pade method 0 , which also does not require a preliminary knowledge of 
the critical index value. One has to differentiate Log of (( 6 }, apply the diagonal Pade 
approximants and define the critical point as the position of the pole nearest to the 
origin. The best estimate obtained this way is xn = 0.919304, with percentage error 
of 1.368%. One can also estimate the value of critical index as a residue 0 , and 
obtain rather disappointing value of 0.73355. 


3.2 Corrected Threshold 

An approach based on the Pade approximants produces the expressions for the 
cross-properties from ”left-to-righf’, extending the series from the dilute regime of 
small x to the high-concentration regime of large x. Alternatively, one can proceed 
from ”right-to-left”, i.e. extending the series from the large x (close to x c ) to small x 

IHlEoJEUi. 

We will first derive an approximation to the high-concentration regime and then 
extrapolate to the less concentrated regime. There is an understanding that physics 
of a 2D high-concentration, regular and irregular composites is related to the so- 
called ’’necks”, certain areas between closely spaced disks mmm. 

Assume also that the initial guess for the threshold value is available from previ¬ 
ous Pade -estimates, and is equal to Xf, = 0.945958. 
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The simplest way to proceed is to look for the solution in the whole region [0,x c ), 
in the form which extends asymptotic expression from l37l , a = a i (x c — x) 1/2 + 
a 2 . This approximation works well for the square lattice of inclusions Cl¬ 
in the case of hexagonal lattice we consider its further extension. 


O = ai(x6~x) S + CC2 +CC3(x6-x) S , ( 10 ) 

where index s is considered as another unknown. All unknowns can be obtained 
from the three starting non-trivial terms of ([ 6 ]), namely (T ~ 1 +2x + 2x 1 + 2x i . Thus 
the parameters equal cq = 2.24674, a 2 = —1.43401, a 3 = 0.0847261, s = 0.832629. 

Let us assume that the true solution a may be found in the same form but with 
exact, yet unknown threshold X c , 

Z = ai(X c -x)~ s + a 2 + a 3 (X c -xy. (11) 

The expression (fTTb may be inverted and X, expressed explicitly, 

^±±y/’ +x . ( 12 ) 

Formula (ITU is a formal expression for the threshold, since E(x) is also un¬ 
known. We can use for E the series in x, so that instead of a true threshold we have 
an effective threshold, X c (x), given in the form of a series in x. For the concrete 
series ©, the following expansion follows. 


= 7 -i/, ( -V(a2-Z) 2 ~4a l a 3 

\ «3 


X r = 2 


X c {x) = x 6 + 0.0134664x 4 + 0.00883052x 5 

+ 0.00647801x 6 - 0.0709217x 7 + 0.0032732x 8 
+ 0.00244442x 9 + 0.00594779x 10 + 0.00482187x* 1 
+ 0.00413887 x 12 +..., (13) 

which should become a true threshold X c as x —> X c . 

Moreover, let us apply re-summation procedure to the expansion ( IT3l ) using the 
diagonal Pade approximants. Finally let us define the sought threshold X* self - 
consistently from the following equations dependent on the approximants order, 

X* n = 0.945958 + 0.0134664x 4 P„,„(X*), (14) 

meaning simply that as we approach the threshold, the RHS of (fl4l> should become 
the threshold. Since the diagonal Pade approximants of the «-th order are defined 
for an even number of terms 2 n, we will also have a sequence of X* n . 

Solving equation ( ITTb . we obtain X* 4 = 0.930222, X* 5 = 0.855009, X* 6 = 
0.9483, X* 7 = 0.932421, X* 8 = 0.946773’z* 9 = 0.941391.±* 10 = 0.94682, X* n = 
0.932752, X* l2 = 0.907423,2f* 13 = 0.903303. The last two estimates for the thresh¬ 
old are good. 
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3.3 Critical index is known 

Also, one can pursue a slightly different strategy, assuming that critical index is 
known (s = 1 /2), and is incorporated into initial approximation. Recalculated pa¬ 
rameters equal ai = 5.12249, U 2 = —5.74972, <23 = 1.52472. For the series ([h}, the 
following expansion follows, 


X c (x) = x 6 — 0.08256 lx 3 + 0.0282108x 4 — 0.000383173x 5 

+ 0.0228241x 6 - 0.0649593x 7 + 0.01561635x 8 
- 0.0091115lx 9 + 0.01874715x 10 + 0.00688507x n 
+ 0.0169516 x 12 +.... (15) 

Let us apply re-summation procedure to the expansion ( fl5l > using super-exponential 
approximants E*(x ) (53]. Finally let us define the sought threshold A* self - consis¬ 
tently, 

X* = 0.945958 -0.08256 lx 3 £* (X* ). (16) 

Since the super-exponential approximants are defined as Ef for arbitrary number 
of terms k, we will also have a sequence of X* k . E.g. 

77 * _ -0.341697* 

E\ — e , 

p* -0.341697e 0 157266i * 
h 2 =e , 

77* -0.341697e 0 157266c5 ' 28382jI + n 

E 3 =e ,..., (1/) 

and so on iteratively. Solving equation (fl6l >. we obtain X* { = 0.901505, X* 2 = 

0.903321, A* 3 = 0.945958, A c * 4 = 0.903404, A c * 5 = 0.916641, X* 6 = 0.903412, 
X* 7 = 0.903556, A* g = 0.903412, X* 9 = 0.903412. 

There is a convergence in the sequence of approximations for the threshold. The 
percentage error achieved for the last point is equal to —0.384537%. 

The method of corrected threshold produces good results based only on the start¬ 
ing twelve terms from the expansion ([6]). in contrast with the Pade-based approxi¬ 
mations, requiring all available terms to gain similar accuracy. The task of extracting 
the threshold, a purely geometrical quantity, from the solution of the physical prob¬ 
lem is not trivial and is relevant to similar attempts to find the threshold for random 
systems from the expressions for some physical quantities El. 

Instead of the super-exponential approximants one can exactly as above apply 
the diagonal Pade approximants, 

A*„ = 0.945958 - 0.082561x 3 P, v ,(A c *). (18) 

Solving equation (fl8l>. we obtain A* 3 = 0.908188, A* 4 = 0.889169, A* 5 = 0.889391, 
X* 6 = 0.887983, A c * 7 = 0.899495, A c * n = 0.903011, X* cV = 0.90296, A* 13 = 
0.9057. 
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Ratio method ©, also works well. It evaluates the threshold through the value 

—+1 

of index and ratio of the series coefficients, x c , n = "<i„ ■ The last point gives rather 

a n— 1 

good estimate, x C: 26 = 0.908801, despite of the oscillations in the dependence on n, 
as seen in Figm 


x c 



Fig. 1 x c calculated by ratio method, compared with the exact threshold. 


4 Critical Index and Amplitude 


Standard way to proceed with critical index calculations when the value of the 
threshold is known can be found in mm. One would first apply the following 
transformation. 


x 

z =- 

X c — X 





(19) 


to the series © in order to make application of the different approximants more 
convenient. 

Then, to such transformed series M\ (z) apply the D — Log transformation and call 
the transformed series M(z). In terms of M(z) one can readily obtain the sequence 
of approximations s n for the critical index ,v, 


s„ = I i m (zPudeA pproximant [M [z}.n.n+ 1]). (20) 

Unfortunately, in the case of © this approach fails. There is no discernible conver¬ 
gence at all within the sequence of s n . Also, even the best result s 12 = 0.573035, is 
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far off the expected 0.5. Failure of the standard approach underscores the need for a 
new methods. 


4.1 Critical Index with D — Log Corrections 

Let us look for a possibility of improving the estimate for the index along the same 
lines as were already employed in the case of a square lattice of inclusions H 1 8H . 
by starting to find a suitable starting approximation for the conductivity and critical 
index. 

Mind that one can derive the expressions for conductivity from ”left-to-right”, 
i.e. extending the series from small x to large x. Alternatively, one can proceed from 
”right-to-left”, i. extending the series from the large x (close to x c ) to small x lfT 8 l| 20 l 
l52l . Let us start with defining reasonable ”right-to-left” zero-approximation, which 
extends the form used in IT8llJ7ll for the square arrays. 

The simplest way to proceed is to look for the solution in the whole region [0,x t ). 
as the formal extension of the expansion, 

o r ~ l = a\(x c -x)~ s + a ,2 + a,T,(x c -x) s + a 4 (x c -x) 2s , (21) 

All parameters in ( fTTT i will be obtained by matching it asymptotically with the trun¬ 
cated series 04 = 1 + 2 x + 2 x 2 + 2 x 3 + 2 x 4 , with the following result, 

°4 '( x ) = (090^^*766 - 5.86967+ 

2.53246(0.9069 - x ) 0 - 520766 - 0.526588(0.9069 -x) 104153 . 

We present below a concrete scheme for calculating both critical index and am¬ 
plitude, based on the idea of corrected approximants m . We will attempt to correct 
the value of ,vo = 0.520766 for the critical index by applying D — Log Pade approx¬ 
imation to the remainder of series ©. 

Let us divide the original series ([ 6 ]) by cr' 4 ~ l (x) given by ( 1221 ). apply to the newly 
found series transformation (fT9l) . then apply D — Log transformation and call the 
transformed series K(z). Finally one can obtain the following sequence of the Pade 
approximations for the corrected critical index, 

s n = so + I i m (z.PadeA pproximant [K[z],n.n + 1]). (23) 

The following ’’corrected” sequence of approximate values for the critical in¬ 
dex can be calculated readily: S 4 = 0.522573, S 5 = 0.518608, s g = 0.554342, 
s 7 = 0.281015, s 8 = —0.209639,, s 9 = 0.279669, si 0 = 0.527055, s n = 0.518543, 
*12 = 0.488502.The last two estimates surround the correct value. 

Generally, one would expect that with adding more terms to the expansion ©, 
quality of estimates for s would improve. As was briefly discussed above, formula 
0 } can be expanded in arbitrary order in x, generating more terms in expansion Q. 
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Fig. 2 Critical index s is calculated by D — Log Corrections method, and compared with the exact 
value. 


Such procedure, of course, is not a rigorous derivation of true expansion, but can be 
used for illustration of the behavior of s n with larger n. 

If y„(z) = PadeA pproximant [K[z].n.n + 1], then 

a n ( x) = ol ~ l M exp ^' l ' Yn (z) dz^j , (24) 

and one can compute numerically corresponding amplitude, 

A n = lim (x c - x) s " a* (x), (25) 

x—>x c 

with Aq = 4.693. Expressions of the type (1241 ) have more general form than sug¬ 
gested before in lfl5l[T6llT9l , based on renormalization methods. 

Convergence for the index above is expected to be supplemented by convergence 
in the sequence of approximate values for critical amplitude, but results are still a 
bit scattered to conclude about the amplitude value. For the last two approximations 
we find An = 4.80599, A 12 = 5.38288, signaling possibility of a larger value than 
4.82, originating from multiplication of the critical amplitude for the square lattice 
by \/3, as suggested by O’Brien l42il . 

To improve the estimates for amplitude A, assume that the value of critical index 
s = 1/2 is given, and construct y„(z) to satisfy the correct value at infinity. There is 
now a good convergence for the amplitude, i.e, in the highest orders, A 10 = 5.09584, 
An = 5.1329, A 12 = 5.14063. Corresponding expression for the approximant 


712 (z) 


bi(z) 

Mz)’ 


( 26 ) 
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where 


bi(z) = —0.079533z 4 — 0.745717z 5 — 2.5712z 6 — 
4.16091z 7 -2.88816z 8 + 0.36028z 9 + 

1.74741z 10 + 0.951728z n — 0.0792987z 12 , 


(27) 


and 


b 2 {z) = 1 + 14.369 lz + 94.745z 2 +380.2z 3 + 
1037.5 lz 4 + 2036.14z 5 + 2961,45z 6 + 

3238. lz 7 + 2667,9z 8 + 1641,88z 9 + 
739.461z I0 + 235.321z n +48.8016z 12 + 
3.81868z 13 . 


(28) 


Corresponding effective conductivity can be obtained numerically. 



(29) 


The maximum error is at x = 0.905 and equals 0.4637%. It turns out that formula 
( l29l ) is good. 


5 Critical Amplitude and Formula for all Concentrations 

For practical applications we suggest below the particular re-summation schemes, 
leading to the analytical expressions for the effective conductivity. 


5.1 Correction with Fade approximates 

Let us ensure the correct critical index already in the starting approximation for 
o' - 1 , so that all parameters in (l30l) are obtained by matching it asymptotically with 
the truncated series 03 = 1 + 2 x + 2 x 2 + 2 x 3 , 



(30) 


3.04972(0.9069 -x) 1 ' 2 - 0.649078(0.9069 - x). 


To extract corrections to the critical amplitude, we divide the original series (] 6 ]» 
by (l30t . apply to the new series transformation ( 1 1 9b . Call the newly found series 
G[z]. Finally build a sequence of the diagonal Pade approximants, so that the ampli¬ 
tudes are expressed by the formula (cq = 5.09924), 


= cq lim (PadeApproximant[G[z],n,n]) 


(31) 


leading to a several reasonable estimates A 7 = 5.26575, An = 5.23882, A 12 = 
5.25781, A 13 = 5.25203. Complete expression for the effective conductivity cor- 
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responding to An can be reconstructed readily. 


<*h{x) = <*l \x)C n (x), 


(32) 


where Ci i(x) = 


ci(x) = 1.15947+ 1.13125 x+1.12212x 2 + 


1.1167X 3 + 3.8727x 4 + 0.824247x 5 - 
2.62954x 6 +1.19135x 7 +1.21923x 8 + 
1.42832x 9 + 1.0608x 10 + 1.53443X 11 ; 


(33) 


and 


c 2 (x) = 1.15947+1.13125 x+1.12212x 2 + 
1.1167x 3 + 3.86892x 4 + 0.849609x 5 - 
2.58112x 6 + 1.11709x 7 + 1.18377x s + 

1.36969x 9 + 1,06062x 10 + x 11 . 


(34) 


Formula (l32l > is practically as good as < [29b . Maximum error is at the point x = 0.905 
and equals 0.563%. 


5.2 Pade approximants. Standard scheme 

Our second suggestion for the conductivity formula valid for all concentrations is 
based on the following conventional considerations Q . Let us first calculate the crit¬ 
ical amplitude A. To this end let us again apply transformation ( IT9l > to the original 
series © to obtain transformed series M\ (z). Then apply to M\ (z) another transfor¬ 
mation to get yet another series, T (z) = M\ (z) -1 ^, in order to get rid of the square- 
root behavior at infinity. In terms of T (z) one can readily obtain the sequence of 
approximations A„ for the critical amplitude A, 



(35) 


There are only few reasonable estimates for the amplitude, Ag = 4.55252, An = 
4.49882, An = 4.64665 and A 13 = 4.68505. The last value is the best if compared 
with the conjectured in If42l . A = 4.82231. 

Following the prescription, the effective conductivity can be easily reconstructed 
in terms of the Pade approximant (corresponding to An) and compared with the 
numerical results in the whole region of concentrations. The maximum error is at 
x = 0.905, and equals —5.67482%. On the other hand, if the conjectured value At, 
is enforced at infinity, through the two-point Pade approximant, the results improve 
and the maximuml error at the same concentration is —3.18511%. Corresponding 
formula for all concentrations, which also respects 24 terms from the series T[z\ is 
given as follows. 
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° p ^ \/0.9069 — x V V 2 (x) ’ 


1.02555 / Vi (x) 


/ 


(36) 


where 


1/, (x) = -0.927562-0.877939x4 0.0406992x 2 + 
0.0440014x 3 + 0.0414973x 4 + 0.0436199x 5 + 

0.319848x 6 4- 0.0110109x 7 - 0.122646x 8 + 
0.0351069x 9 + 0.0439523x 10 + 0.0380654x n + 
1.01499 x 12 +x 13 


(37) 


and 


V 2 (x) = -1.07571 +2.09854x-2.17 187 x 2 + 
2.23064x 3 -2.3122x 4 + 2.374X 5 - 
2.1397x 6 + 1,87791x 7 - 1.78516x 8 + 

1,86446x 9 - 1,94838x 10 + 2.03264X 1 r - 

x 12 


(38) 


Various expressions are shown in Fig{3] Note, that significant deviations of the 
Corrected Pade formula (1441 and of the Standard Pade formula ( 1361 ) from the refer¬ 
ence rational expression ©. start around x = 0.85. All formulas start to depart from 
the original series around x = 0.8. The two formulae, (l44l > and ( l36l >. happen to be 
very close to each other almost everywhere, except in the immediate vicinity of the 
critical point. 


5.3 Accurate final formula 

According to our calculations, based on various re-summation techniques applied 
to the series ©, we conclude that the critical amplitude is in the interval from 5.14 
to 5.24, by 6 — 9% higher than following naively to O’Brien’s 4.82. 

Below we present an exceptionally accurate and more compact formula for the 
effective conductivity ( l32l > valid for all concentrations. 

Let us start with modified expression ( l30l > taking into account also the O’Brien 
suggestion already in the starting approximation for the amplitude in o' All re¬ 
maining parameters in ( l30l > are obtained by matching it asymptotically with the 
truncated series o 2 = 1 + 2 x + 2 x 2 , 



(39) 


2.13365(0.9069 -x ) 1 / 2 -0.328432(0.9069 -x). 


Repeating the procedure developed in subsection 15.11 we receive several rea¬ 
sonable estimates for the critical amplitude, A 7 = 5.18112, An = 5.15534, A 12 = 
5.19509,Ai3 = 5.18766. 

Complete expression for the effective conductivity corresponding to the first es¬ 
timate for the amplitude, is given as follows 
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a^x) = G r 2 - l {x)F 1 {x), (40) 

and Fi(x) = 444, where 

v ' AM 

/i(x) =52.0141 + 10.3198x-38.8957x 2 + 4.70555x 3 + 

4.89777x 4 + 4.6887x 5 + 0.476241x 6 + 7.49464x 7 , ( } 


and 

fi (x) = 52.0141 + 10.3198x- 38.8957 x 2 + 2.17078x 3 + 

5.80088x 4 + 6.03946x 5 + 1.80866x 6 + x 7 . ( ’ 

The formulae predict a sharp increase from o 2 (0.906) = 166.708, to (0.9068) = 
513.352, in the immediate vicinity of the threshold, where other approaches 1421 l4l 
[9), fail to to produce an estimate. At the largest concentration x = 0.9068993 men¬ 
tioned in (42), the conductivity is very large, 8375.34. This formula (|40| ) after slight 
modifications can be written in the form ([5]). 

Asymptotic expression can be extracted from for the approximant (l40t . 


ct 


* 


5.18112 

V0.9069-X 


-6.229231. 


(43) 


Even closer agreement with numerical results of m is achieved with approxi¬ 
mant corresponding to A 13 . 


CTi* 3 (x) = 02 r \x)F n {x), 


(44) 


where F u (x) = 

/i(x) = 1.49313+ 1.30576 x + 0.383574x 2 +0.4677 13x 3 + 

0.471121x 4 + 0.510435x 5 + 0.256682x e + 

0.434917x 7 + 0.813868 x 8 + 0.961464x 9 + C 

0.317194x 10 + 0.377055x u - 1.2022 x 12 -0.931575x 13 ; 

and 

/ 3 (x) = 1.49313 +1.30576 x + 0.383574x 2 + 0.394949x 3 + 

0.44785x 4 + 0.503394x 5 + 0.303285x 6 + 

0.271498x 7 + 0.732764x 8 + 0.827239x 9 + (46) 

0.25509x 10 + 0.239752X 11 - 1,26489x 12 - x 13 . 

It describes even more accurately than ( l40l >. the numerical data in the interval 
from x = 0.85 up to the critical point. The maximum error for the formula (l44l) is 
truly negligible, —0.042%. 

Asymptotic expression can be extracted from for the approximant (l44t . 


a 


* 


5.18766 
V0.9069 — x 


6.2371. 


(47) 
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cr 



Fig. 3 Our formula d441 (solid) is compared with the standard Pade approximant d361 (dotted) and 
rational approximation 0 (dashed). The series 0 is shown with dashed line. 


5.3.1 Role of randomness 

For random two-dimensional composite we obtained recently ll36l . the following 
closed-form expression for the effective conductivity. 


a*(x) = 0.121708/^ r (x)x 

(0.64454-r—1.38151 ).v+0.72278 


exp 


(jc—0.9069)- 


jt(j+0.435329)+0.3582 
(x-0.9069) 2 


- 0.815613 sinh 


(- 


0171(*+0.494058)'\ 
x—0.9069 


(48) 


with 

^ f , (0.419645x4- l) 3 - 45214 

/o ’ r(X) “ VI-1.10266* 


(49) 


Closed-form expression for the effective conductivity of the regular hexagonal 
array of disks is given by (l44l >. Since the two expressions are defined in the same 
domain of concentrations, a comparison can explicitly quantify the role of a ran¬ 
domness (irregularity) of the composite. In order to estimate an enhancement factor 
due to randomness we use ratio of (l48l > to (144b . In particular, the enhancement factor 
at x = 0.906, is equal to 104.593. In Figure[4] such an enhancement factor is shown 
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R(x) 



Fig. 4 Ratio R(x) = ^ ^ of the effective conductivity for the random composite to the effective 
conductivity of the hexagonal regular lattice calculated with <48b and ED, respectively. 


in the region of high concentrations. 


6 Interpolation with High-concentration Pade approximants 


When two expansions © and ( 143[ i are available, the problem of reconstruction 
greatly simplifies and can be solved upfront in terms of Pade approximants. 

This approach requires as an input at least two parameters from weak and 
strong- coupling (high-concentration) regimes, including the value of amplitude 
A = 5.18112 from (03). Similar problem for random composites was considered 
in H. 

Assume that the next-order term, B = —6.22923 from (03), is known in advance. 
The high-concentration limit. In terms of ^-variable ( fl9) . the strong-coupling limit 
is simply 

(T~ -^Vz + B + 0{z~ l/2 ). (50) 

\ x c 


The Pade approximants all conditioned to give a constant value as z —> 0 are given 
below. 


P2,l(z) 

P3,2(z) 

P4,3(Z) 


M 1+ ^ + ^) 

1 +fr;7j+lr+ftz _3/2 +§) 

\ + p 5 J^ + £i.+fc z -y2 


(51) 


The unknowns in (f5Tb will be obtained by the asymptotic conditioning to (150) and 
©. In all orders j3 = -^=. Explicitly, in original variables, the following expressions 
transpire. 
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P2,l(x) 


/- X -i 4.9348 

V 0.9069-a- ' 0.9069-v” 


■4.11284 


V O.0O69-.V + 1 -3 2 856 


_ (0-608173v / roro^+ 1 - 26S63 >+°- 677749 v / 0.‘)0^--v + 1 13282 

P3 ’ 2[X> ~ -(0- 747 325 v ^=;+l)r+0.677749 v 4^+1.13282 


(52) 


5.4414^ 1+3.76414— 7 = 4 - 7 ^@i+1.9'7396z(j:) - 3 / 2 +24Z68ll 

/ \ _ _ v V _ \AM ; W _ z(*) z 

/4 ’ 3( ' 3 k+4!90893^ 

The approximants are strictly non-negative and respect the structure of ©, e.g. 
for small x, 

Pa^{x) ~ 1 +2x+2x z + 0(x i ), (53) 

since all lower-order powers generated by square-roots, are suppressed by design. 
But in higher order, emerging integer powers of roots should be suppressed again 
and again, to make sure that only integer powers of x are present. As x —> x c , 


P4 3 (x)~A(x c -x) 1/2 + fi + 0((x c -x) 1/2 ), 


(54) 


and only integer powers of a square-root appear in higher-orders. Both P 3 . 2 ( x ) and 
Pa.} (x) give good estimates for the conductivity, from below and above respectively. 
Their simple arithmetic average works better than each of the approximants. The 
bounds hold till the very core of the high-concentration regime, till x = 0.906. 

Particularly clear form is achieved for the resistivity, an inverse of conductivity, 
r(z) = (p{z))~\ e.g„ 


r 3,4 (z) = 


3.76815 +1.97396^ +0.902145z + 0.183776z 3 / 2 
3.76815 + 1.97396 yjz + 7.7368 k + 3.76414z 3 / 2 + z 2 


(55) 


With the variable X = yfz, the resistivity problem is reduced to studying the 
sequence of Pade approximants R n = r n n t i (A), n = 1,2..i/2, with X G [0,°°), and 
analogy with the Stieltjes truncated moment problem | I][l4][2), is complete as 
long as the resistivity expands at X —> °° in the Laurent polynomial with the sign- 
alternating coefficients, coinciding with the “Stieltjes-moments” /i/ ; (see e.g., 

[48l . were the original work of Stieltjes is explained very clearly). 

The moments formally define corresponding Stieltjes integral as X —> °°, 


i: 


’ d(j)(u) 
u -\- X 


E(- i)V*x - * -1 - 


-o(x~ l ), 


(56) 


k=0 


/ is even m, and /Ik = f^Li k d(j)(u). Approximant R„{X ) should match l l56l > asymp¬ 
totically. 

The Stieltjes moment problem can possess a unique solution or multiple solu¬ 
tions, dependent on the behavior of the moments, in contrast with the problem of 
moments for the finite interval 101251 05), which is solved uniquely if the solution 
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exists l50l . The role of variable is played by the contrast parameter, while in our 
case of a high-contrast composite, the variable is X. 

In our setup,there are just two moments available and resistivity is reconstructed 
using also a finite number of coefficients in the expansion at small X. I.e., the re¬ 
duced (truncated) two point Pade approximation is considered, also tightly related 
to the moment problem 041(291 l27l . In fact, even pure interpolation problem can 
be presented as a moment problem. We obtain here upper and lower bounds for 
resistivity (conductivity) in a good agreement with simulations t42l . 

It does seem interesting and non-trivial that the effective resistivity (conductivity) 
can be presented in the form of a Stieltjes integral I50ll49ll48l . when the variable 
(O is used. 


6.1 Independent estimation of the amplitude B 

We intend to calculate the amplitude B independent on previous estimates. Start 
with the choice of the simplest approximant as zero-approximation, 




(58) 


The way how we proceeded above was to look for multiplicative corrections to 
some plausible ’’zero-order” approximate solution. We can also look for an additive 
corrections in a similar fashion. To this end subtract (l58l ) from (0 to get some new 
series g{x). Change the variable x =y 2 to bring the series to a standard form. The 
diagonal Pade approximants to the series g(y) are supposed to give a correction to 
the value of 1, suggested by (l58l >. To calculate the correction one has to find the 
value of the corresponding approximant as y —> fx~ c . The following sequence of 
approximations for the amplitude B can be calculated now readily. 


B n = 1 + PadeApproximant[g(y —>■ y/x^),n,n\. 


(59) 


The sequence of approximations is shown in FigJ3 
??? move B in the y-axis 

There is clear saturation of the results for larger n, and Bo6 = —5.94966. One can 
reconstruct the expression for conductivity corresponding to B 2 6 in additive form 


°26 (x)=P 1,0 (x) + F 26 ( x ) , 


(60) 


where F 26 (x) = 
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B 



Fig. 5 Sequence of approximations B„ calculated from 1591 . 


F 2 (x) = -5.71388^- 

1,5564x - 0.358877x 3 / 2 - 2.18519x 2 + 0.0918426x 5 / 2 - 1,59468x 3 - 
0.149418x 7 / 2 — 1.47691x 4 — 0.366848x 9 / 2 — 1.49733x 5 -0.56432x n / 2 - 
1.58738x 6 + 0.21344x 13 / 2 — 1.3108 lx 7 —0. 366156a 15 / 2 + 15.1037x 8 — 
15.1703x 17 / 2 -6.38227x 9 + 0.576004x 19 / 2 -2.5147x 10 + 0.715526x 21 / 2 - 
1.53752X 11 + 0.28655x 23 / 2 - 1.1985 lx 12 + 5.951 lx 25 / 2 + 0.55801 lx 13 , 

and 

F 6 (x) = 1 + 0.622415 v^x — 0.27066x + 0.294568x 3 / 2 — 0.00182913x 2 + 
0.0875453x 5 / 2 + 0.0832152x 3 + 0.098912x 7 / 2 + 0.114562x 4 + 0.116471x 9 / 2 + 

0.133685x 5 + 0.133737X 11 / 2 - 0.0205003x 6 + 0.045100lx 13 / 2 + 

0.134146x 7 - 2.77482x 15 / 2 + 1.6976x 8 + 3.12806x 17 / 2 -0.87267x 9 + 
0.16541x 19 / 2 -0.179645x 10 + 0.0404152x 21 / 2 + 0.000620289x n + 

0.0514796x 23 / 2 -0.986857 a 12 -0.58415x 25 / 2 + 0.388415x 13 . 

(62) 

The maximum error for the formula (l60t is very small, 0.0824%, only slightly in¬ 
ferior compared with (l44l >. The amplitude B is firmly in the interval (5.95,6.22), 
according to our best two formulae. 


7 Discussion of the ansatz (121130b 


In the case of a square lattice of inclusions li42l |30l [3T1 5] [18], we looked for the 
solution in a simple form. 


cr[ z = ai(x c -x) 1 / 2 + a 2 , 


c 


7Z 

4’ 




(63) 
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and obtained the unknowns from the two starting terms of the corresponding series, 
C7~ l+2x + 2x 2 + 2x 3 + 2x 4 + ... (64) 

3/2 

Then, CX\ = « 2.784, (X 2 = (1 — 7l), same form as obtained asymptotically in 

03 , with exactly the same value for the leading amplitude as obtained in |[32l . 

Also the lower bound for amplitude A is equal to 2.753, and the upper bound is 
equal to = 2.856, as can be found directly from the corresponding corrected Pade 
sequences for the critical amplitude m. 

Formula (l63l > despite its asymptotic nature, turned out to be rather accurate in the 
whole region of concentrations. We try to understand below why it is so. 

Let us subtract the approximant (l63l ) from the series ( l64l >. apply to the new series 
transformation ( flbl ). Then we apply to such transformed series another procedure, 
intended to find corrections to the values of amplitudes a 1 and CD- Such task is non¬ 
trivial, especially when one is interested in analytical solutions. It can be solved 
using general form of root approximants derived in Il20ll52l . 

Oadd = boz 2 {{b\z+ 1) S1 +b 2 z 2 ) 52 (65) 

under asymptotic condition 

Oadd -d l ^z + d 2 , asz-^-oo. ( 66 ) 

Elementary power-counting gives si —3/2, s 2 = —3 /4. All other unknowns can 
now be determined uniquely in a standard fashion from the condition of asymptotic 
equivalence as z —> 0. Final expression 


®add — 


0.0556033x 2 


10 7 R^QS - A 2 ( j-69302* 2 + < 1-98243.V , i \ 3 / 2 A 3/4 ’ 

(U./833yo X) y (0.785398 —x)- + V0.785398 -x + / 


(67) 


can be re-expanded in the vicinity of x c with the result 
0.0184973 


6 add — 


V0.785398 -x 


-0.0118315 +0(y/x^x), 


( 68 ) 


indicating only small corrections to the values of amplitudes. Such asymptotic sta¬ 
bility of all amplitudes additionally justifies the ansatz, and final corrected expres¬ 
sion a sq = O’, 1 + O a dd, appears to be just slightly larger than ( 1 6 31 . Note that mod¬ 
ified Pade approximants as described above, are only able to produce additive cor¬ 
rections in the form o a dd — d^z + 0(-^) as z —> °°- 

In the case of hexagonal lattice, such simple proposition as |63] does not ap¬ 
pear to be stable in the sense described above. We have to try lengthier expres¬ 
sions of the same type, such as (l30l >. Additive correction in the form (7,,,/d = 
boz 4 (b 2 z 2 + (b\z + l) 3 / 2 ) ~ 7 / 4 , or 
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Gadd ( t ) 


0.0022082 lx 4 




leads to the very small, almost negligible asymptotic corrections to the ansatz 
(l30l) . E.g., the leading amplitude changes to the value of 5.09925. Such asymp¬ 
totic stability of all amplitudes justifies the ansatz. Of course, it also appears to be 
reasonable when compared with the whole body of numerical results. The lower 
bound= 5.0925, and the upper bound bound= 5.298, can be found directly from the 
corresponding corrected Pade sequences. 


8 Square and Hexagonal United 

From the physical standpoint of there is no qualitative difference between the prop¬ 
erties of hexagonal and square lattice arrangements of inclusions. Therefore one 
might expect that a single expression exists for the effective conductivity of the two 
cases. 

Mathematically one is confronted with the following problem: for the functions 
of two variables <J sq (x,x s c q ) and ai iex {x 1 x h c ex ), to find the transformation or relation 
which connects the two functions. (Here x% ex = x c ). 

Assuming that the expressions for both lattices are different only with respect 
to lattice parameters simplifies the task, but is not necessary. The problem is really 
simplified due to similar leading asymptotic terms in the dilute and highly concen¬ 
trated limits. On general grounds, one can expect that up to some simply behaving 
’’correcting” function of a properly chosen non-dimensional concentration, the two 
functions are identical. Below we do not solve the problem from the first principles, 
but address it within the limits of some accurate approximate approach. 

We intend to express O sq and 0] lex in terms of the corresponding non-dimensional 

x q —x x! wx —x 

variables, Z sq = -S— and Z/, CT = c , lex , respectively. Each of the variables is in the 

X c X! c 

range between 0 and 1. 

Then, we formulate a new ansatz which turns to be good both for square and 
hexagonal lattices. 



where u is a control parameter introduced by the self-similar renormalization lf22l 
|24l applied to the asymptotic form (l63l >. One can obtain the unknowns from the 
three starting terms of the corresponding series, which happen to be identical for 
both lattices under investigation. 

Then, the method of (15.11) . when the ansatz (f70l> is corrected through application 
of the Pade approximants, is applied. Emerging diagonal Pade-sequences for criti¬ 
cal amplitudes are convergent for both lattices and good results are simultaneously 
achieved in the same order, employing 24 terms from the corresponding expansions. 
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We select from the emerging sequences only approximants which are also holo- 
morphic functions. Not all approximants generated by the procedure are holomor- 
phic. The holomorphy of diagonal Pade approximants in a given domain implies 
their uniform convergence inside this domain (A.A. Gonchar, see l26l ). 

Corresponding corrective Pade approximants, Cor^*, Cor sq 2 , are given below, 
compatible with the formulae given below. For the hexagonal lattice, 

oUZ) = o c - hex {Z)Cor l [?{Z), (71) 

and for the square lattice, 

a; q {Z) = a c ^(Z)Cor\l(Z). (72) 

The initial approximation for the hexagonal lattice, 

, , 4.5509 (1 - 0.637832 y/Z) 149198 

a c,he* ( Z ) = -i- - Z-J. -, (73) 

\JZ 

and for the square lattice. 


ct c «(Z) 


3.29343 (1 - 0.659155^) 11074 

Vz 


Correction term for the hexagonal lattice. 


Corft(Z) 


and for the square lattice, 


cor\ ex {Z ) 
cor'^ x (Z) ’ 


Cor\l(Z) 


cor\ q (Z) 
cor 2 q (Z )' 


(74) 


(75) 


(76) 


Numerators and denominators of these expressions are given by polynomials, 

cor'{ ex (Z) = 4.77682 - 2.70811Z-20.9607Z 2 + 139.454Z 3 - 

411.652Z 4 + 752.321Z 5 - 935.979Z 6 + 824.735Z 7 - 520.883Z S + (77) 

232.778Z 9 - 70.2807Z 10 + 12.8747Z 11 - 1.07872Z 12 


cor\ ex (Z) =4.07267 + 1.71405Z-35.8279Z 2 + 172.528Z 3 — 

46L417Z 4 + 801.75Z 5 - 964.958Z 6 + 828.769Z 7 - (78) 

512.137Z 8 + 224.545Z 9 —66.7017Z 10 + 12.0594Z 11 —Z 12 

cor q (Z) = 10.0658 - 55.5091Z+96.8247Z 2 - 149.99Z 3 + 

227.446Z 4 - 298.272Z 5 + 346.9Z 6 - 341.722Z 7 + 258.688Z 8 - 
138.017Z 9 +47.9727Z 10 — 9.69297Z 11 +0.860976Z 12 


( 79 ) 
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cor«(Z) = 10.3657 - 57.833Z+ 105.045Z 2 - 175.546Z 3 + 

298~.234Z 4 —438.138Z 5 + 532.347Z 6 — 508.883Z 7 + (80) 

362.562Z 8 — 182.102Z 9 + 60.161Z 10 — 11.6563Z 11 +Z 12 

The ratio of final expressions for the conductivity of corresponding lattices, 
^SrplrT- ■ can be plotted (as Z] xex = Z sa = Z), as shown in Fig. [6] 

°sq\ Zj sq) * 


he X 



Fig. 6 The ratio of final expressions for the conductivity of corresponding lattices, > can 

be plotted (as Z) wx = Z sq = Z). 


It turns out that the ratio is bounded function of Z, and changes monotonously 
from 1 (Z = 1) to 1.669 (Z = 0). The last number is not far from the O’Brien sug¬ 


gestion ll42l . and is simply 
critical amplitudes. 



= 5.08318, = 2.834284, are the 


9 Dirichlet summation. Large-rz behavior of series coefficients. 

We will try to evaluate how the coefficients of the series behave at large-n. From 
the practical viewpoint it is beneficial to have such information (if available), to be 
included into resummation procedure. The so-called Borel summation is known to 
render filed-theoretical calculations more consistent. With a similar goal, we employ 
the ordinary Dirichlet’s series, defined conventionally 0(c) = JT" , a n n ~ c , where a n 
stands for the coefficients of the original series. 

Essential difference distinguishes the general theory of Dirichlet’s series from 
the simpler theory of power series. The region of convergence of a power series is 
determined by the position of the nearest singular points of the function which it 
represents. The circle of convergence extends up to the nearest singular point. No 
such simple relation holds in the general case of Dirichlet’s series. When convergent 
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in a portion of the plane they only may represent a function regular all over the plane, 
or in a wider region of it. 

However in an important case relevant to our study, the line of convergence nec¬ 
essarily contains at least one singularity. It is covered by the following theorem: 

Theorem 10 lf28l . If all the coefficients of the series are positive or zero, then the 
real point of the line of convergence is a singular point of the function represented 
by the series. 

We conjecture, following |:35l, that for large-/; the sum-function of coefficients, 
S n = ci\+ci 2 +... +a n , behave as follows. 


S n ~ 8n c ' log £ (zz). 

Then, Dirichlet’s series can be written explicitly in the form lf35l . 

<j) (c) = 8cr (e + 1) (c - ci) 1 + g (c), 


(81) 


(82) 


where g(c) stands for the regular part, and 8 is a parameter. This expression is valid 
at c > ci, where the Dirichlet’s series are convergent. 

In order to return to the physical region of variables x and conductivity, let us 
apply the following transformation 


with the inverse 


c{x) 

x(c) 


x c (x + x c ) 
x c -x 

x c c - x c 2 
x c +c 


(83) 

(84) 


with ci = x c . 

The singular part of the conductivity after such transformation is expressed in the 
form 


2-‘->sr( £ +i)(,+*)fe) 

CT s (x) = ---4 


(85) 


and we should also set e = —1/2. Parameter 8 is simply connected with the critical 
amplitude A, 8 = —A —. 

r \/2nx c 

Finally, 



This expression should also be regularized at small x, so that 


(5 S ,r{x) 



( 87 ) 


Close to critical point it can be expanded. 
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Os,r{x) — 


y/x c ~: 


A A{— X + X c ) 

2^~ Wfi +0(fe-») / ). 


( 88 ) 


After extracting the singular part from the series, the regular part expands for 
small x into the following expression (only few low-order terms are shown) 


g{x) ~ 1 - 


3A^/x 

4x c 


+ 2x — ■ 


7 Ax 3 / 2 
16x c 2 


+ 2 x 2 + 0(x 5/2 ), 


(89) 


which is an expansion in ^fx. To this expansion we apply the diagonal Pade approx- 
imants. Presence of fractional powers can be easily taken into account by change 
of variables, x = y 2 , leading to doubling the number of approximants which can be 
constructed, compared with series of only integer powers. 

E.g., in the lowest orders, in addition to a standard polynomial ratio with inte- 
ger highest power, _ 0 . 460544;c+0 . 3 91 4 1 4 V ^ + p there is another ratio a4667 ^ +1 , with 
fractional highest power, which can be considered as a diagonal Pade approximant 
too. Only the former-type polynomial ratios will be presented below, since the latter- 
type ratios do not bring better results in the current context. 

Our goal now is to calculate the second, constant term in expansion close to x c , 
denoted above as B. The correction to the constant term in the expansion emerges 
directly from the Pade approximant calculated at x = x c , 


Bn — 


2y/x~c 


+ PadeApproximant [g(x —> x c ),n,n]. 


(90) 


We receive several reasonable estimates for the amplitude B: B$ = —6.40157, 
B 6 = —6.28506,^7 = —6.27028,B 8 = -6.33762 ,B n = -6.29595, B n = -6.29695, 
= -6.29842. 

Explicitly in 7-th order, 


of* = G sr (x) + PadeA pproximant [g [x], 7,7]. 
Corresponding expressions for the singular part of solution, 

3 A) 


Oi,r(x) = 


2 ^ 33 / 4 ^ 


(91) 


( 92 ) 


and for the regular part g(x) = given by the Pade approximant, we find 

Gi (x) = 23.7835 - 88.5524^ - 39,6443x + 71,3743x3/ 2 + 
36.2957x 2 - 12.3254x 5 / 2 -5.54733x 3 - 1.28303x 7 / 2 - 
4.81208x 4 - 1.00508x 9 / 2 -4.36713x 5 - 1.3826X 11 / 2 — 

6,04028x 6 + 4.92363X 13 / 2 + 0.518137x 7 ; 
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G 2 (x) = 23.7835 + 13.3695^- 29.9175x- 18.0154x 3 / 2 + 

8.21262x 2 + 6.49253x 5 / 2 + 0.387922x 3 + 1 ,2539 x 7 / 2 + 

0.50647x 4 + 0.689113x 9 / 2 + 0.505114x 5 + 0.554584x‘ V 2 - 
1.31685x 6 —0.498839 x 13 / 2 + x 7 . 

The maximum error for the formula is small, just —0.1602%. 

The formulae predict the following values, af*(0.906) = 166.494, af’(0.9068) = 
512.7472, O® (0.9068993) = 8376.58. These value are very close to the predictions 
already presented above. 

We conclude that our conjecture concerning the large-« behavior of the sum- 
function of the coefficients, is in a good agreement with available numerical data. 
Also the estimates for B, which stem from the conjecture, is close to other estimates 
from the present paper. 

Algorithms and mathematical methods used above, are based on asymptotic 
power-series for the effective conductivity and various resummation techniques to 
ensure their convergence. Such approach is typical for Computational Science of 
Composite Materials. It is interesting to compare such approach and classic method¬ 
ology based on direct solutions of PDE’s. 


10 Application of Lubrication Theory 

To find the effective conductivity in a classic way, one has to consider the local 
problem for Laplace equation describing regular hexagonal lattice of cylindrical 
inclusions. Such a study can be based on the Lubrication theory [.121, applicable for 
an asymptotic regime of large, ideally conducting inclusions. It has to be applied 
in conjunction with some averaging technique to derive effective conductivity. It is 
expedient first to consider inclusions with finite conductivity A, and then to consider 
the limit A —> °°. 

Main idea of the Lubrication theory consists in replacing the original bound¬ 
ary problem with another, corresponding to a simpler geometry (see LigJT}. I.e. the 
original hexagonal elementary cell is replaced by a circle off radii b. Using so-called 
’’fast” variables (|, rj) and the corresponding local polar coordinates (r, 0) we arrive 
at the following problem (for details see GDI ED) 



(95) 



dr dr 

u = 0, r = b, 


(96) 


(97) 


where a is the radius of inclusions. Lor definiteness, the external flux is taken in 
such a way that the macroscopic flow is presented by the potential liq = (xi ,x 2 ) and 
the flux by the vector (1,1) (for details see |30l ). The problem (|95|>-(|97 |i has the 
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1 



Fig. 7 The hexagonal cell with the disk of the radius a is approximated by the circle cell of the 
radius b. 


solution 


Ahrcos 0 + N 2 r sin 0, r<a , 

[M\V cos 0 -T {M 2 v -f sin 0, a<r<b, 


(98) 


where the constants are determined by the boundary conditions 

- At, - (A-l)(fc 2 -fl 2 ' 


N, =n 2 = 


[6 2 +a 2 -A(6 2 -a 2 )] ’ 


A/i = Mi = — 


(A-l)c/ 2 


=^2 = - 


[f> 2 +a 2 -A(6 2 -a 2 )]’ 
(A—1)A 2 


(99) 


According to the Lubrication approach fl2l l5ll, let us consider an external contour 
for the cell, as a circle of varying radii 





( 100 ) 


Integration is conducted over the quarter of the elementary cell, shown in Fig[8] 
Following general prescriptions of the averaging method, we derive averaged coef¬ 
ficient 


a = 


1 

W\ 



du 



du 



, ( 101 ) 


where |f21 = 2\/3- The integration is performed to satisfy also the relation (1 100b : in 
particular, /?(<!;) is considered as a corresponding functions of varying radius. 
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Fig. 8 Approximation of the hexagonal cell by the circle cell of the variable radius b(%). 


10.1 Lubrication Approximation and Correction 

After some transformations we receive the following expression for the effective 
conductivity (or thermal conductivity) as the function of the inclusion size a. 


(2\/3a 2 ) tan^ 1 ( —r==^ 
_/ \ v ' V3Vl-a 2 

o(a) = - 


'\-a l 


+ 1 + 


i(y/3a 2 ) (f-|sin 1 (^)) 


. 4V3a 2 
3-v/1 -a 2 


tan 


5 tan -1 


/3a y/fa 2 1j \j\—a 


a +1 y 

2(-V / 3fl 2 ^T+v / 3a-1) V^l-" 2 


I tarn 


-v/3a 2 -l(V3fl+a-2)+(l+V3)a(l-V3a)+2 

1 ( 2y/T^ 


, , / 3a 2 +2./3 (3a 2 —1) +2 

- 4—3a 2 


( 102 ) 


As the inclusion size tend to its limiting value, a —>• 1, the leading term in the con¬ 
ductivity of the ideally conducting inclusions can be found in the familiar form. 


do ^ 


y/l—a 


(103) 
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When expressed in terms of volume fraction of inclusions (11031 ). coincides with 
Keller’s formula ((3}. 

The first (constant) correction term to the formulae ( 1 1 02b . can be also obtained, 
leading to ’’shifted” expression for the conductivity in the critical region, 

cti ~ C7 0 -5.10217. (104) 

Formula (11021 ) works rather well, with accuracy less than 2%, for concentrations as 
low as x ~ 0.82. Its predictions for concentrations very close to x c , are also very near 
to predictions from other formulae given above (see Figj9]). Formula (1102b becomes 
invalid forx < 0.3023. 


cr 



Fig. 9 a calculated by formulas §5} (solid line), by {47} (dotted line) and by fl04l (dashed line). 


In the case of a square lattice of inclusions. Lubrication theory gives the follow¬ 
ing asymptotic result da, 

3/2 

CT W-^= -1 - (105) 

Formula f 105 1 should be compared with the more accurate result of I37l . ll38l . 

7T 3 /2 

a(x) " Vf^~ 7r+L (106) 

It appears that Lubrication theory assumptions, concerning reduction of the el¬ 
ementary cell to a circle, work better for the hexagonal lattice than for the square 
lattice. In both cases the correction term is overestimated. 

Classic approach to PDE’s solution thus is limited to high-concentration asymp¬ 
totic regime with strong interactions between inclusions. 
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On the other side, the whole well-developed family of self consistent methods 
which include Maxwell’s approach, effective medium approximations, differential 
schemes etc., are valid only for a dilute composites when interactions between in¬ 
clusions do not matter mu. 

In contrast, computational methods of the present paper are applicable every¬ 
where. 

Let us derive an interpolation formula by matching the two limiting expres¬ 
sions. © and (1 104b . The method of sewing the two limiting behaviors together 
will be chosen to employ the main idea of Section 0 First we assume that the 
high-concentration formula (1 1 04b holds everywhere and then derive an additive cor¬ 
rection in the form of the diagonal Pade approximants in such a way that also the 
low-concentration limit ([ 6 ]) is respected. It turns out that such approach not only gen¬ 
erates another good interpolation formula, but also calculates an additive correction 
to the amplitude B. Technically, one should only replace the expression (l 86 l > with 
(1 1 04b and extract it from the [ 6 l leading to the new series g(x) and to corresponding 
approximations to the sought amplitude, 

B n = —5.10217 + PadeApproximant[g(x ^ x c ),n,n]. (107) 

We receive several reasonable estimates for the amplitude B: B$ = —6.37811, 
B 6 = —6.29179,B 7 = —6.28019, B & = -6.42952, B n = -6.29702,B 12 = -6.29908, 
B 13 = —6.32249. These results are only slightly higher than estimates obtained 
above in Section [9] Interpolation formula corresponding to B-j is as accurate as its 
counterpart suggested in Section[9] 


11 Random composite from hexagonal representative cell 

In the present paper, the numerical computations for random composites are per¬ 
formed for the hexagonal representative cell. The number of inclusions per cell can 
be taken arbitrary large, hence the shape of the cell does somewhat influences the 
final result. 

The hexagonal lattice serves as the domain Q, where random composite is gener¬ 
ated as a probabilistic distribution of disks of radius r (particles), by means of some 
Monte-Carlo algorithm (protocol) El- 

Algorithm 1, random sequential addition (RSA). First random point is randomly 
distributed in Q. Second point is randomly distributed in Q with exception of the 
small circular region of radius r surrounding the first point. Hence, the distribution 
of the second random point is conditional and depends on the first random point. 
More points, up to some number N, can be generated, conditioned that circular re¬ 
gions around all previous points are excluded from Q. This joint random variable 
for all points correctly determines sought probabilistic distribution. But the com¬ 
puter simulations work only up to concentrations as high as 0.5773, hence is the 
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main RSA limitation. To overcome the limitation and to penetrate the region of 
larger concentrations, one has to apply some extrapolation technique. 

Algorithm 2, random walks (RW) employed also in |36| . N-random points are 
generated, at first being put onto the nodes of the hexagonal array. Let each point 
move in a randomly chosen direction with some step. Thus each center obtains new 
complex coordinate. This move is repeated many times, without particles overlap. 
If particle does overlap with some previously generated, it remains blocked at this 
step. After a large number of walks the obtained locations of the centers can be 
considered as a sought statistical realization, defining random composite. 

RW protocol can be applied for arbitrary concentrations including those very 
close to x c , which stands also for the maximum volume fraction of random compos¬ 
ites. At x = xc = -=, we arrive to the regular hexagonal array of disks. 

The effective conductivity of random composite is also expected to tend to infin¬ 
ity as a power-law, as the concentration x tends to the maximal value x c , 

<t(x) — A(x c —x) — L (108) 

The superconductivity critical exponent s believed to be close to | ~ 1.3 l46ll . 
much different from the regular case. The critical amplitude A is an unknown non- 
universal parameter. We demonstrate below that s depends on the protocol, and sug¬ 
gest simple way to decrease the dependence on protocol. Still, more studies are 
needed with different protocols. 

Algorithm 2 allows to obtain the following power series in concentration, 

0 RW = 1 + 2x + 2x 2 + 4.2372 lx 3 + 6.8975x 4 . (109) 

The higher order polynomial representations fail to give a non-zero value for the 
fourth-order coefficient. 

Reasonable estimate for the critical index s can be obtained already from the 
D — Log formula combined with the transformation ( 1 1 9| i. (l20l) . Namely, the result is 
S 2 = 1.43811, and for the amplitude we obtainA = 1.21973. 

The algorithm 1 produced the following series in concentration l36l . 

a RSA = i + 2x + 2x 2 + 5.00392x 3 + 6.3495x 4 . (110) 

The coefficients on x^ (k = 5.6.7,8) vanish in (11 101) with the precision 10 l0 . 

Good estimate for the critical index s can be obtained already from the D — Log 
formula ( l20l >. The results are S 2 = 1.28522 for the critical index, and A = 1.57678 
for the amplitude. 

Ideally, we would like to have s and A to be evaluated independent on protocol, 
but can hope only that combining two different protocols can decrease the depen¬ 
dence of s on protocols, because errors of the two protocols can compensate. 

Assume that both schemes should lead to the same index, amplitude and thresh¬ 
old. Let us form a simple product. 


a J = \J a RW a RSA . 


(Ill) 



Title Suppressed Due to Excessive Length 


33 


11.1 D-Log estimates 


Apply now the D — Log technique combined with the transformation (fl9l ). to the 
series (11111) . The result is s = S 2 = 1.34715, better than for each of the individual 
components. 

Simple addition of (11 101) and (11091 ) also leads to a good estimate 1.34888, by the 
D — Log technique. 

Slightly better result is achieved for the geometrical mean of the series, 

j'jRW-RSA 

m _ zo ° 

° a RW + a RSA’ 


and s = S 2 = 1.34542. The coefficients in the expansion for small x, 

1+2x + 2a: 2 + 4.62056x 3 + 6.6235x 4 , (113) 

are formed as a compromise between the two algorithms. 

The effective conductivity can be reconstructed mum no, from an effective 
critical index (or -function). After some calculations, we obtain 


°f to = 

( 0.9069 
V 0.9069—a 


3.24319e°' 441389tan ‘ 
-0.5l5l66) im ° 9 ( 


(2.18756 + ^S) X 

x(x+0.0245056)+0.176696 
(0.9069 —a) 2 


0.00466513 


(114) 


Also, the critical amplitude evaluates as 1.423. Eq. (11141 ) works as good as any other 
formula for the effective conductivity obtained in |[36l . 


11.2 ’’Single pole” approximation 


Critical index can be estimated also from a standard representation for the derivative 


B a (x) = d x \og(o M (x)) ~ —— 
x c -x 


(115) 


as x —> x c , thus defining critical index as the residue in the corresponding single 
pole. 

Outside of the immediate vicinity of the critical point a diagonal Pade approxi- 
mant is assumed for the residue estimationj6l, but such approach fails in the case 
under study. Let us use another representation, in the form of a factor approximant 

m, 


B a (x) 


2 {p2X T 1) S2 


(116) 


with the following values for parameters b 2 = 7.84091, so = —0.140629, found for 
the series (II 12b . 

Formula (11161 ) leads to the simple expression for the critical index 
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S = 2x c (b 2 X c + \) S2 , 


(117) 


and to the value s = 1.35129. The effective conductivity can be reconstructed as 
follows, 

C7*(x) = 

2K^b 2 x+iy2+i 2Fl (i’ s 2 + i;^ + 2;ifc^)) 

(nb 2 +2V3)(s 2 +l) 

through the hypergeometric function. The ’’single pole” approximation (1115b is in 
fact equivalent to the particular case of the hypergeometric function. 

For the RSA-series (11 1 Ob . the same approach gives s = 1.31786, while for the 
RW-series (1 109b . s = 1.37978. The difference between the two algorithms is small 
compared to all others methods employed for the index estimations. 


(118) 



11.3 Corrected Index. Scheme 1 


We follow below the general idea of |T8]|36), also explained in Subsection 14.11 At 
first, one should obtain an approximate solution explicitly as a factor approximant 
1231 El. Then we attempt to correct the form of the initial approximation with 
additional factor, originated from the part of series which did not participate in the 
formation of the initial approximation, following literally the way leading to (f24t . 

The simplest factor approximant can be calculated, 

(x -)- i ) 1-04882 

/oW _ (1 - 1 . 10266 jt ) 0 - 862622 (119) 

Such approximant satisfy the two non-trivial starting terms from the series (11 131) . 
and incorporates the accepted value of the threshold x c . It predicts for the critical 
index the value .vq = 0.862622. 

In the next step we attempt to correct sq using the D — Log-correction approach 
mm, as described also in Subsection 14. II Let us form the following ratio, -fij-- 
Repeating the same steps that lead to the corrected expression for the index [23] 
we obtain the corrected value sj = 1.32067. Corresponding amplitude is equal to 
1.48267. 

The conductivity can be reconstructed in a closed form. Calculating correspond¬ 
ing integral with B -function II 1 911T61 -f*2,3 (z), 




5.71085z 2 

12.4679z 3 + 10.3351z 2 + 4.31945z+ 1’ 


( 120 ) 


we obtain 
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C7*(x) = 1.77719 


(*+i)‘ 


(1 — 1.10266a) 0 - 862622 


g 0.465101 tan-* 


„ t0.545059.t+0.412586t 0444818 i 
* V 0.9069 —x I ' 


+0.0241854.V+1) 1807- 

(0.9069—jc) 2 , 


( 121 ) 


11.4 Corrected Index. Scheme 2 


Let us start from the initial approximation (111 9b . and recast it more generally as 

/oW = (1 - (122) 

where R(x) stands for the regular part of i ll 191 ). In what follows we attempt to correct 
/o(x) differently than above, assuming instead of sq some functional dependence 
S(x). 

As x —1 x c , S(x) —i s c , the corrected value. The function S(x) will be designed in 
such a way, that it smoothly interpolates between the initial value ,vo and the sought 
value s c . The corrected functional form for the conductivity is now 

fW = (l-i)-%(r). (123) 

X c 

From (1 1 23b one can express S(x'), but only formally since f* (x) is not known. But we 
can use its asymptotic form ( 111 3b . express S(x) as a series and apply some resumma¬ 
tion procedure (e.g. Pade technique). Finally calculate the limit of the approximants 

as x —>• x c . 

In what follows the ratio C(x) = R \ > , stands for an asymptotic form of the 
singular part of the solution, and as x —> 0 


S{x) 


log(C(x)) 

'og(l--j)' 


(124) 


which can be easily expanded in powers x, around the value of so- It appears that 
one can construct a single meaningful Pade approximant, 

4.91072x 2 + 0.703479x + 0.862622 
X} ~ 3.00966x 2 + 0.815512x+1 ’ 


and find the corrected index, s c = S(x c ) = 1.31426. Now we also possess a complete 
expression for conductivity (1 1 23b . 

Scheme 2 due to its simplicity, can always lead to the analytical expression. But 
Scheme 1 seems to be the most flexible. It also turns out to be weakly dependent 
on the starting approximation fo(x). Indeed, if another starting approximation is 
considered, 

(2x+l) 0 - 355391 

/oW = (1 — 1.10266x) 116919 ’ (126) 
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the corrected index remains good, S 3 = 1.31094. 

The conductivity again can be reconstructed in a closed form. Calculating corre¬ 
sponding integral with /3-function ^ 34 ( 2 ), 


P3A(z) 


9.85652z 3 +4.06592z 2 

69.5337Z 4 + 35.9326z 3 + 20.6483z 2 + 6.17673z + 1 ’ 


(127) 


we obtain rather lengthy expression, 

( 2 *+l ) 0 - 355391 


o 3 *(x) = 1.05615 

exp (0.0615685tan- 1 (2.23114 - -0.0789727tan" 1 (4.693 - 0 ^ 358 t )) 

,, /r(l-0.245701)+0.138583\ 0081481 , /x(j:+ 0.415099)+0.099469 n \ 00106051 
x V (0.9069 — v) 2 ) ' V 


(0.9069-jc) 2 


J 


(128) 

The form of expressions (1121b . (1 1 28b is unlikely to be guessed as an independent 
approximant. 


12 Conclusion 

Based on estimates for the critical amplitudes A and B, we derived an accurate and 
relatively compact formula for the effective conductivity (l40t valid for all concen¬ 
trations, including the most interesting regime of very high concentrations. For the 
high-concentration limit, in addition to the amplitude value of 5.18112, we deduce 
also that the next order (constant) term B, equals —6.22923. It is possible to ex¬ 
tract more coefficients in the high-concentration expansion based on the formula 
®. Dirichlet summation is suggested to extract an arbitrary large-/; behavior of 
the coefficients. 

When two expansions around different points ([ 6 | and ( 1431 are available, the prob¬ 
lem of reconstruction can be solved in terms of high-concentration Pade approxi- 
mants, implying that the effective resistivity (conductivity) can be presented in the 
form of a Stieltjes integral, in terms of the variable X=. —. Such Pade approx- 

\f X C X 

imants give tight lower and upper bounds for the conductivity, valid up to the very 
high x. 

Such properties as the superconductivity critical index and threshold for con¬ 
ductivity, can be calculated from the series (] 6 ]». In the case of truncated series, the 
standard Pade approximants are not able to describe the correct asymptotic behav¬ 
ior in the high-concentration limit, where in addition to the leading critical exponent 
also a non-trivial sub-leading exponent(s) plays the role fl20l[52 1. On the other hand 
when such a non-trivial asymptotic behavior is treated separately with different type 
of approximants, the Pade approximants are able to account for the correction. Such 
patchwork approximations appear to be more accurate and powerful than approxi¬ 
mating conventionally with a single type of approximants. 













Title Suppressed Due to Excessive Length 


37 


Simple functional relation between the effective conductivity of the hexagonal 
and square lattices is suggested, expressed in terms of some bounded monotonous 
function of a non-dimensional concentration of inclusions. Getting an accurate for¬ 
mula in this case, means that correct asymptotic behavior (l43l > is indeed can be 
extracted from the series ©, and together they determine the behavior in the whole 
interval with good accuracy. Neglecting the high-concentration regime dominated 
by necks, is not admissible. 

We also considered a classic approach based on Lubrication theory and con¬ 
cluded that it can be applied strictly within the high-concentration asymptotic 
regime. In contrary, the celebrated Maxwell’s approach, effective medium approxi¬ 
mations and differential schemes are valid only for a dilute composites ED- Com¬ 
putational approach and results of the present paper are applicable everywhere. 

We conclude that approach based on the long power series for the effective con¬ 
ductivity as a function of particle volume fraction can be consistently applied in the 
important case of highly conducting (superconducting) inclusions. Based on our in¬ 
vestigation we put forward the final formula ©, for the effective conductivity of the 
hexagonal array. 

Acknowledgements The authors are grateful to Leonid Berlyand for stimulating discussion. 
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